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Abstract 

Apoptosis is a highly regulated cell death mechanism involved in many physiological processes. A key 
component of extrinsically activated apoptosis is the death receptor Fas, which, on binding to its cognate 
ligand FasL, oligomerize to form the death-inducing signaling complex. Motivated by recent experimental 
data, we propose a mathematical model of death ligand-receptor dynamics where FasL acts as a clus- 
tering agent for Fas, which form locally stable signaling platforms through proximity-induced receptor 
interactions. Significantly, the model exhibits hysteresis, providing an upstream mechanism for bistability 
and robustness. At low receptor concentrations, the bistability is contingent on the trimerism of FasL. 
Moreover, irreversible bistability, representing a committed cell death decision, emerges at high concen- 
trations, which may be achieved through receptor pre-association or localization onto membrane lipid 
rafts. Thus, our model provides a novel theory for these observed biological phenomena within the unified 
context of bistability. Importantly, as Fas interactions initiate the extrinsic apoptotic pathway, our model 
also suggests a mechanism by which cells may function as bistable life/death switches independently of 
any such dynamics in their downstream components. Our results highlight the role of death receptors in 
deciding cell fate and add to the signal processing capabilities attributed to receptor clustering. 

Author Summary 

Many prominent diseases, most notably cancer, arise from an imbalance between the rates of cell growth 
and death in the body. This is often due to mutations that disrupt a cell death program called apoptosis. 
Here, we focus on the extrinsic pathway of apoptotic activation, which is initiated upon detection of an 
external death signal, encoded by a death ligand, by its corresponding death receptor. Through the tools 
of mathematical analysis, we find that a novel model of death ligand-receptor interactions based on recent 
experimental data possesses the capacity for bistability. Consequently, the model supports threshold- 
like switching between unambiguous life and death states — intuitively, the defining characteristic of an 
effective cell death mechanism. We thus highlight the role of death receptors, the first component along 
the apoptotic pathway, in deciding cell fate. Furthermore, the model suggests an explanation for various 
biologically observed phenomena, including the trimeric character of the death ligand and the tendency 
for death receptors to colocalize, in terms of bistability. Our work hence informs the molecular basis 
of the apoptotic point-of-no-return, and may infiuence future drug therapies against cancer and other 
diseases. 

Introduction 

Apoptosis is a coordinated cell death program employed by multicellular organisms that plays a central 
role in many physiological processes. Normal function of apoptosis is critical for development, tissue 
homeostasis, cell termination, and immune response, and its disruption is associated with pathological 
conditions such as developmental defects, neurodegenerative disorders, autoimmune disorders, and tu- 
morigenesis [lji5j. Due to its biological significance, much effort has been devoted to uncovering the 
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pathways governing apoptosis. Indeed, recent progress has enabled the prohferation of mathematical 
models, both mechanistic and integrative [e.g., IBHH], which together have offered profound insights into 
the underlying molecular interactions. The current work takes a similarly mathematical approach and 
hence inherits from this legacy. 

There are two main pathways of apoptotic activation: the extrinsic (receptor-mediated) pathway and 
the intrinsic (mitochondrial) pathway, both of which are highly regulated [151 Ull • this study, we focus 
on the core machinery of the extrinsic pathway, which is initiated upon detection of an extracellular death 
signal, e.g., FasL, a homotrimeric ligand that binds to its cognate transmembrane death receptor. Fas 
(CD95/Apo-l), in a 1:3 ratio. This clusters the intracellular receptor death domains and promotes the 
ligation of FADD, forming the death-inducing signaling complex (DISC) pTHTO] . The DISC catalyzes 
the activation of initiator caspases, e.g., caspase-8, through death effector domain interactions. Initiator 
caspases then activate effector caspases, e.g., caspase-3, which ultimately execute cell death by direct 
cleavage of cellular targets [2QH23] . 

Apoptosis is typically viewed as a bistable system, with a sharp all-or-none switch between attract- 
ing life and death states. This bistability is important for conferring robustness [24]. Consequently, 
researchers have used computational models to identify and study potential sources of bistability in 
apoptosis, including positive caspase feedback [8 , inhibition of DISC by cFLIP [2, cooperativity in 
apoptosome formation [10 , double-negative caspase feedback through XIAP [IT], and double-negative 
feedback in Bcl-2 protein interactions [25 . In this work, we propose that bistability may be induced 
upstream by the death receptors themselves. 

The current model of death ligand-receptor dynamics assumes that FasL activates Fas by direct 
crosslinking, producing a DISC concentration that varies smoothly with the ligand input [26 . However, 
recent structural data [27] suggests a different view. In particular. Fas was found in both closed and 
open forms, only the latter of which allowed FADD binding and hence transduction of the apoptotic 
signal. Moreover, open Fas were observed to pair-stabilize through stem helix interactions. This affords a 
mechanism for bistability, similar to the Ising model in ferromagnetism [28], where open Fas, presumably 
disfavored relative to their native closed forms [29], are able to sustain their conformations even after 
removal of the initial stimulus promoting receptor opening, past a certain critical density of open Fas. 
This induces hysteresis in the concentration of active, signaling receptors and therefore in apoptosis. 

We studied this proposed mechanism by formulating and analyzing a mathematical model. The 
essential interpretation is that FasL acts as a clustering platform for Fas, which establish contacts with 
other Fas through pairwise and higher-order interactions to form units capable of hysteresis (Figure [T]). 
At low receptor concentrations, the model exhibits bistability provided that the number of receptors that 
each ligand can coordinate is at least three. This hence gives a theory for the trimeric character of FasL. 
Furthermore, at high concentrations, for example, through receptor pre-association [3QH32] or localization 
onto lipid rafts [33], irreversible bistability is achieved, implementing a permanent cell death decision. 
Thus, our model suggests a primary role for death receptors in deciding cell fate. Moreover, our results 
offer novel functional interpretations of ligand trimerism and receptor pre-association and localization 
within the unified context of bistability. 

Results 

Model formulation 

Constructing a mathematical model of Fas dynamics is not entirely straightforward as receptors can 
form highly oligomeric clusters [271 ES]. A standard dynamical systems description would therefore 
require an exponentially large number of state variables to account for all combinatorial configurations. 
To circumvent this, we considered the problem at the level of individual clusters. Each cluster can be 
represented by a tuple denoting the numbers of its molecular constituents, the cluster association being 
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implicit, so only these molecule numbers need be tracked. 

In our model, a cluster is indexed by a tuple (L, X, Y, Z), where L represents FasL and X, y, and Z 
are three posited forms of Fas, denoting closed, open and unstable, and open and stable, i.e., active and 
signaling, receptors, respectively. Within a cluster, we assumed a complete interaction graph and defined 
the reactions 
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The first reaction describes spontaneous receptor opening and closing; the second, constitutive destabi- 
lization of open Fas; the third, ligand-independent receptor cluster-stabilization; and the fourth, ligand- 
dependent receptor cluster-stabilization (Figure [2| . The orders of the cluster-stabilization events are 
limited by the parameters m and n, which capture the effects of receptor density and Fas coordination 
by FasL, respectively. Although only pair-stabilization (m = n = 2) has been observed experimentally 
[27] , higher-order analogues, for example, as facilitated by globular interactions, are not unreasonable. 

Formally, these reactions are to be interpreted as state transitions on the space of cluster tuples. 
However, the reaction notation is suggestive, highlighting the contribution of each elementary event, which 
we modeled using constant reaction rates (for simplicity, we set uniform rate constants /cs*^ and for all 
ligand-independent and -dependent cluster-stabilization reactions of molecularity respectively). Then 
on making a continuum approximation, we reinterpreted the molecule numbers as local concentrations 
and applied the law of mass action to produce a dynamical system for each cluster in the concentrations 
(/, X, z) of (I/, X, F, Z). Validity of the model requires that the molecular concentrations are not too low 
and that the timescale of receptor conformational change is short compared to that of cluster dissociation. 

To study the long-term behavior of the model, we solved the system at steady state (denoted by the 
subscript oo). Introducing the nondimensionalizations 
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where s is a characteristic concentration and t is time, and 

(3a) 
(3b) 

z = 2, . . . , m, (3c) 



ko 
kc 

kc 



nl^) = \±, i = 2,...,n, (3d) 



this is 



^oo /^o^oo, (4b) 

where a = -\- r] -\- ( is the nondimensional total receptor density, and ("oo is given by considering 

f = E E ^'c- E ^ + ^ E E ^'c-^ E - (5) 

i=2 j=l /c=l i=2 j = l k=l 

and solving d^/dr = with (^, ?7, C) ^ (^oo ? ^oo ? Coo ) ? a polynomial in Coo of degree max{m,n}. Clearly, 
the model is bistable only if max{m, n} > 3 (two stable nodes must be separated by an unstable node as 
the model is effectively one-dimensional in Q. 

We used C as a measure of the apoptotic activation of a cluster. In principle, all open receptors 
contribute to apoptotic signaling, but r] is small, at least at steady state (since <C 1 due to the 
assumed prevalence of the closed form [29]), and so can be neglected. 



Bistability and receptor clustering 

While n measures the coordination capacity of FasL and hence may be equated with its oligomeric 
order (e.g., n = 3 in the biological context), an appropriate value for m, relating to the total receptor 
concentration, is somewhat more elusive. Therefore, we began our analysis by performing a simple 
receptor density estimate. Approximating the cell as a cube of linear dimension ^ 10 /im, the associated 
volume of ~ 1 pL implies the correspondence 1 nM ^ 600 molecules ~ 10~^ molecules/nm^ on restricting 
to the membrane, i.e., by averaging over the surface area of ~ 600 /im^. Thus, for a conservative receptor 
concentration estimate of 100 nM [71 [9l |T2l , the number of Fas molecules in the neighborhood of each 
receptor is only ~ 1, assuming a charateristic size of 100 nm. We hence found that receptors may be very 
sparsely distributed. In this low density mode, high-order Fas interactions in the absence of ligand can 
be neglected (m = 2). Therefore, in this context, bistability is possible only if n > 3, and the trimerism 
of FasL thus demonstrates the lowest-order complexity required for bistability. 

From the form of d(/dr^ this bistability is reversible as a function of the FasL concentration A since 
the governing polynomial for Coo is of degree only m = 2 at A = 0. This suggests that at the cluster 
level, the cell death decision can be reversed, which may have adverse effects on cellular and genomic 
integrity. However, irreversible bistability at higher receptor densities may also be achieved. Researchers 
have observed tendencies for death receptors both to pre-associate as dimers or trimers [3QH32] and to 
selectively localize onto membrane lipid rafts [33^ . The result of either of these processes may be to increase 
the local receptor concentration. In this high density mode, we set m > 3, as the preceeding approximation 
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is no longer valid. Irreversible bistability then becomes attainable, representing a committed cell death 
decision. 

For the remainder of the study, we incorporated both the low and high receptor density regimes into a 
single model by setting m = 3, using a as a continuous transition parameter. Furthermore, we set n = 3 
to correspond to observed biology. 

Characterization of the steady-state surface 

Calculation of the steady-state activation curves showed that the model indeed exhibits bistability (Figure 
[3| for reasonable parameter choices (Methods). Thus, we established the possiblity of a novel bistability 
mechanism in extrinsic apoptosis. The associated hysteresis enables threshold switching between well- 
separated low and high activation states. Biologically, these define local signals of life and death, which 
are integrated at the cell level to compute the overall apoptotic response. 

As per the previous analysis, reversibility of the bistability is dependent on a, with irreversibility 
emerging for a sufficiently high. This suggests a bivariate parameterization of Coo, producing a multivalued 
steady-state surface over (A, <j)-space (Figure [4|. The result is a cusp catastrophe, an elementary object 
of catastrophe theory, which studies how small perturbations in certain parameters can lead to large and 
sudden changes in the behavior of a nonlinear system [34 . A more instructive view of the dependence of 
the model's qualitative structure on A and a is shown in Figure |5j 

Sensitivity and robustness analyses 

We then focused on the activation and deactivation thresholds A±, respectively, defining the bistable 
regime. These are the points at which the steady state switches discontinuously from one branch to 
the other, and are given by the values of A at which the hysteresis curve turns, i.e., at dX/d^oo = 
(Figure |6|. We performed a sensitivity analysis of X± by measuring the effects of perturbing the 
model parameters about baseline values (Methods). For each threshold-parameter pair, we computed a 
normalized sensitivity S by linear regression. 

Strong effects of a, and were observed (Figure [t]); for the corresponding Fas thresholds (± 

(2) — [3) 

at A = A±, respectively, the parameters a, a^o, k^i , and were emphasized. Thus, the bistability 
thresholds do not appear particularly robust. However, the data reveal that essentially all parameter 
sets sampled were bistable. This suggests a weaker form of robustness, namely, robustness of bistability, 
which nevertheless supports life and death decisions over a wide operating range. 

To probe this further, we sampled parameters with increasing spread D about baseline values and 
computed the fraction / of parameter sets that remained bistable (Methods). The results show that / 
has an exponential form (Figure [8|. Extrapolating to I) ^ oo, the data suggest an asymptotic bistable 
fraction of foo ^ 0.4. Hence, robustness of bistability remains substantial even under significant parameter 
variation. 

Cell-level cluster integration 

Thus far, we have considered only the apoptotic activation of an individual cluster. To obtain the more 
biologically relevant cell-level activation, we must integrate over all clusters. In principle, this integration 
should account for intercluster transport as well as any intrinsic differences between clusters, e.g., as due 
to spatial inhomogeneities. Here, however, we provide as demonstration only a very simple integration 
scheme. Specifically, we assumed that clusters are identical (apart from their parameter values, which 
are drawn randomly) and independent, and that FasL is homogeneous over the cell membrane. Then we 
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can express the normalized cell activation as 

C"(A) = ^^|^, (6) 

where the subscript i denotes reference to cluster i. 

A characteristic cell- level hysteresis curve is shown in Figure |9] As is immediately evident, such 
integration is a smoothing operator, averaging over the sharp thresholds of each cluster. Thus, the 
cell-level signal may be graded even though its constituents are not. Note, however, that the lack of a 
sudden switch from low to high Fas signaling does not necessarily imply the same at the level of the 
caspases which ultimately govern cell death, as downstream components may possess switching behaviors 

laiiiioiiiiiES]. 



Model discrimination 



Finally, we sought to outline protocols to experimentally discriminate our model against the prevail- 
ing crosslinking model [26 , which we considered in a slightly simplified form [35 . To be precise, the 
crosslinking model that we used has the reactions 
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where L is FasL, R is Fas, and Ci is the complex FasLiFas^ for i = 1,2,3. With 
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(continuing the notational convention that lowercase letters denote the concentrations of their uppercase 
counterparts), the steady-state solution under mass-action dynamics is 
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and 

A = A + 7i + 72 + 73, (11a) 
cr = p + 7i + 272 + 373 (lib) 

are the total ligand and receptor concentrations, respectively. In analogy with our proposed model, 
hereafter called the cluster model, we used 

C = 71 + 272 + 373 = a - p (12) 

as a measure of the apoptotic signal. 



Hyperactive mutants 



Clearly, the crosslinking model has only one steady state, while the cluster model is capable of bistability. 
This hence provides a ready discrimination criterion. Although tracing out the associated hysteresis 
curve may be problematic, we can nevertheless probe for bistability by using hyperactive mutants, e.g., 
the mutation of He 313 to Asp in Fas, which stabilizes the open conformation and enhances apoptotic 
activity [27 . 

Specifically, we considered an experimental setup in which the concentrations of FasL and Fas, both 
wildtype and mutant, can be controlled, and in which the apoptotic signal can be measured, e.g., through 
the degree of FADD binding or of caspase activation. Hence we can map out the response curves at 
various levels of mutant penetrance. Denoting mutant Fas by Za (we assumed that mutant Fas cannot 
close, so there is no distinction between the stable and unstable open forms), we characterized the mutant 
penetrance by the mutant population fraction A = Ca/^, where Ca = ^a/^ is the nondimensional mutant 
Fas concentration and a = a + Ca is the total receptor concentration, composed of contributions from 
both wildtype (a) and mutant (Ca) forms. We assumed no other functional differences between wildtype 
and mutant Fas. 

Proceeding first for the crosslinking model, at fixed a, the amount of Fas bound by FasL is determined 
only by A. Hence we assumed that a fraction (p = (p{A) of all receptors is bound. Since wildtype and 
mutant Fas are functionally identical by assumption, the fraction bound for each of the wildtype and 
mutant populations is also (p{A) by independence of the recruitment process. Therefore, under the 
crosslinking model, varying A at fixed a yields an invariant response curve for the active wildtype Fas 
fraction (p = (^/a. 

In contrast, for the cluster model, we expected mutant receptor cluster-interactions to affect the 



wildtype response. Accordingly, the reactions (Ic) and (Id) were amended for interaction with Za by 
replacing with 



jY -^kZ^{i-j-k)ZA^ U -k')Y^{k^ k') Z^{i-j-k) Za, 



l(0 



L^jY^kZ^{i-j-k)ZA^L^ {j - /cO r + (/c + k') Z^{i-j-k) Za, 
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respectively. This gives the analogue 

f = E E E ^'cr-' E fc' + ^ E -i' E E c'^a--'^ E fc' - ««c (m) 

2=2 j = l k=0 k' = l 2=2 j = l /c=0 k' = l 

of (§. 

As seen in Figure [lOj receptor interactions indeed cause the apoptotic signal to increase with A even 
after accounting for the effect of mutants. This is because mutants can activate wildtype receptors by 
pushing the cluster past its switching threshold. Furthermore, the convergence to the active cluster state 
at high A provides evidence for bistability. Thus, the variance of the (^-response curve at various A can 
be used for model discrimination. 



Steady-state invariants 

Alternatively, if working with mutants should prove difficult, we provide also a discrimination test based 
on steady-state invariants, i.e., functions that vanish at steady state. Clearly, for each model, ( = d(/dr 
provides a steady-state invariant since C = necessarily at steady state. However, the difficulty lies 
in expressing ( solely in terms of variables that are experimentally accessible. For example, current 
technology may not allow the concentrations 7^ to be measured accurately, if at all. Therefore, all such 
variables must be eliminated. Rate constants were considered parameters and so were not subject to this 
rule. 

We assumed the same experimental setup as above and hence expressed each model invariant in terms 
of A, cr, and ("oo, giving functions of the form ("(A, a, ("oo; ci), where a encompasses all model parameters 
(Methods). The task then is to use |C|, with (A, a, Coo) provided by experiment, to assess the fit of a 
model. However, the parameters a remain unknown, so this assessment cannot proceed directly. Instead, 
we considered the best possible fit min^ |C| over all parameters. A high value of min^ |C| indicates a poor 
best-case fit and hence that a model is unlikely to be correct. Clearly, prior knowledge of a can be used 
to guide the invariant to biologically plausible fits. 

To demonstrate that model discrimination using steady-state invariants is practical, we generated 
synthetic data from each model, calculating the accessible concentrations (A, a, Coo) ^ot each parameter 
set. This gives two sets of model-generated data. For each data set, we computed the best-fit invariant 
error e = min^ rms \ (\ for each model, where rms is the root mean square operator. The results suggest 



that this test can correctly identify the model from the data (Figure 11). 

The systems that we have presently considered are simple enough that experimentally inaccessible 
variables can be eliminated by hand. For more complicated systems, the tools of computational algebraic 
geometry, notably Grobner bases, may prove useful; for such an application, see [36^ . 



Discussion 

In this work, we showed through analysis of a mathematical model that receptor clustering can support 
bistability and hysteresis in apoptosis through a higher-order analogue of biologically observed Fas pair- 
stabilization [27 . Hence we add to the signal processing activities in which receptor clustering has been 
suggested to participate [37-39 . This bistability plays an important functional role by enabling robust 
threshold switching between life and death states. Significantly, our results indicate potential key roles 
for ligand trimerism [TT and receptor pre-association [30^^3? and localization onto membrane lipid rafts 
[33] . Thus, we provide novel interpretations for these phenomena within the unified context of bistability. 

Our model suggests an additional cell death decision, supplementing those that have been studied 
previously [71 [SI [TOl [HI [25] . Critically, the proposed decision is implemented upstream at the very death 
receptors that initially detect the death signal encoded by FasL. This decision is therefore apical in that it 
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precedes all others in the system. Consequently, it operates independently of all intracellular components 
and so offers a general mechanism for bistability, even in cell lines with, for example, only feedforward 
caspase-activation networks [7, [131 [14 . Thus, receptor cluster-activation may explain how an effective 
apoptotic decision is implemented in such cells. Moreover, this suggests a novel target for induced cell 
termination in the treatment of disease [T. 

We believe that our model provides an attractive theory for the observed biology. Although unlikely 
to be correct in mechanistic detail, the model may nevertheless reflect reality at a qualitative level. 
The significance of our work hence lies in its capacity to guide future research. We therefore readily 
invite experiment, which can reveal the true nature of the molecular mechanisms involved. Given their 
structural and functional homology, similar investigations on other members of the tumor necrosis factor 
receptor family may also prove fruitful. Such work serves to further our understanding of the formation 
and mode of action of complex signaling platforms such as the DISC, which in this view may be considered 
the macromolecular aggregates of active Fas. 

Methods 
Parameter selection 

The rationale for the choices m = n = 3 is presented in the text; here, we further defend these by noting 
that no new behaviors are introduced with m or n > 3. The remaining parameter values were guided by 
the following considerations. Specifically, we required hCq and /^^^ <C 1 due to the assumed stabilities of 
the receptor species; all other parameters were assumed to be close to 0{1). Within these constraints, 
parameters were selected to ensure that A+ is of the correct order of magnitude [71 [3 [121 [13]. The 
baseline parameter values used were a = 1, /t:o = 2 x 10~^, hiu = 10~^, /^s^^ = 0.1, i^i^"^ = 0.5, = 1, 
and /^p'' = 5. 

Parameter sampling 

To analyze the effects of variability in the model parameters, parameter values were sampled from a 
log-normal distribution, characterized by a variation coefficient defined as the ratio of the standard 
deviation to the median of the distribution. For the sensitivity analysis, 500 parameter sets were drawn 
at D = 0.25; for the robustness analysis, 250 parameter sets were drawn over < D < 5; and for the 
cell-level integration, 100 parameter sets were drawn at D = 0.25. All parameters were drawn about 
baseline median values. 

Sensitivity analysis 

For each threshold-parameter pair, linear regression was performed on the threshold data against the 
parameter data, each normalized by reference values. For parameters, the reference is the baseline 
(median) value; for thresholds, the reference is the threshold computed at baseline parameters. The 
normalized sensitivity S was defined as the slope of the linear regression. 

Steady-state invariants 

The cluster invariant was derived by considering ( at steady state, i.e., with (^, ^ (^oo,^oo), and 
identifying A A. Similarly, the crosslinking invariant was derived by considering ( at steady state and 
identifying poo ^ cr — Coo- For the full forms of the invariants, see Protocol SI. 

For the model discrimination computation, 100 parameter sets (A, a) were drawn from log-normal 
distributions with median (0.2, 1) at D = (0.2,0.5). The active Fas concentration was calculated for 
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each parameter set for each model at basehne parameters {k: = 0.1 for the crosshnking model); for the 
cluster model, if bistability was observed, one of the stable values of ("oo was chosen at random. The 
invariant error e was minimized using SLSQP with a lower bound of 10~^ for all parameters. 

Computational platform 

All calculations were performed with Sage 4.5 [40^, using NumPy/SciPy [IT for numerical computation 
and matplotlib [42] for data visualization. The Sage worksheet containing all computations is provided in 
the Supporting Information (Protocol SI) and can also be downloaded from http://www.sagenb.org/ 
|home/pub/1224/| or |http: //www. courant .nyu. edu/~ho/[ 
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Figure 1. Cartoon of model interactions. The transmembrane death receptor Fas natively adopts 
a closed conformation, but can open to allow the binding of FADD, an adaptor molecule that facilitates 
apoptotic signal transduction. Open Fas can self-stabilize via stem helix and globular interactions, 
which is enhanced by receptor clustering through association with the ligand FasL. 
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Figure 2. Schematic of cluster-stabilization reactions. Examples of ligand-independent 
cluster-stabilization reactions involving unstable {Y) and stable (Z) open receptors of molecularities 
two (A), three (B), and four (C). Higher-order reactions follow the same pattern. Ligand-dependent 
reactions are identical except that FasL (L) must be added to each reacting state. 
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Figure 3. Steady-state activation curves. The steady-state active Fas concentration shows 
bistability and hysteresis as a function of the FasL concentration A (stable, solid lines; unstable, dashed 
lines). At low receptor concentrations cr, the bistability is reversible, but irreversibility emerges for a 
sufficiently high, representing a committed cell death decision. All parameters set at baseline values 
unless otherwise noted. 
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Figure 4. Steady-state activation surface. The steady-state surface for the active Fas 
concentration ("oo as a function of the FasL and total Fas concentrations A and a, respectively, is folded, 
indicating the existence of singularities, across which the system's steady-state behavior switches 
between monost ability and bistability (stable, blue; unstable, red). All parameters set at baseline values 
unless otherwise noted. 
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Figure 5. Steady state diagram. Steady state diagram identifying the regions of parameter space 
supporting monostability (colored) or bistability (gray) as a function of the FasL and total Fas 
concentrations A and a, respectively. The monostable region is colored as a heat map corresponding to 
the steady-state active Fas concentration Coo- Irreversible bistability is indicated by the extension of the 
bistable region to the axis A = 0. 



Bistability in Apoptosis by Receptor Clustering 



15 




FasL (A) 



Figure 6. Bistability thresholds. The activation (red) and deactivation (blue) thresholds X± 
characterizing the bistable regime (green) are defined as the concentrations A of FasL at which the 
steady-state active Fas concentration (oo (black) switches discontinuously from one branch to the other 
(stable, solid line; unstable, dashed line). 
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Figure 7. Sensitivity analysis of bistability thresholds. The robustness of the bistability 
thresholds is investigated by measuring the effects of perturbating the model parameters about baseline 
values. For each threshold-parameter pair, a normalized sensitivity S is computed by linear regression. 
Top, sensitivities for the FasL thresholds A±; bottom, sensitivities for the corresponding Fas thresholds 
C± at FasL concentrations A = A±, respectively. 
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Figure 8. Robustness of bistability. The fraction / of parameter sets that exhibit bistabihty as a 
function of the samphng variabihty D follows the exponential form / = /oo + (1 — foo)^~^^^^^ where 
/oo is the asymptotic bistable fraction. The fitted value of /oo ^ 0.4 suggests that this robustness 
remains substantial even as I) ^ oo. 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



FasL (A) 




Deactivation threslnold (A ) 



Figure 9. Cell-level cluster integration. The apoptotic signals of all Fas clusters are integrated to 
produce a normalized cell activation < C^^^ < 1. The resulting hysteresis curve on Co©^^ ^ function 
of the FasL concentration A is graded due to the heterogeneity of the bistability thresholds \± across 
the clusters (top). Despite this variability a strong linear dependence persists between \± (bottom; the 
valid region > A_ is shown in green). 



Bistability in Apoptosis by Receptor Clustering 



17 




0.00 



0.05 



0.10 



0.15 0.20 
FasL (A) 



0.25 



Figure 10. Model discrimination using hyperactive mutants. The wildtype response curve, 
giving the steady-state active wildtype Fas fraction as a function of the FasL concentration A (stable, 
solid lines; unstable, dashed lines), of the cluster model varies with the mutant population fraction A, 
reflecting receptor interactions absent in the crosslinking model. The total receptor concentration is 
fixed at a = 1. All parameters set at baseline values unless otherwise noted. 
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Figure 11. Model discrimination using steady-state invariants. Steady-state invariants are fit 
to synthetic data generated from each model. For each model-data pair, the invariant error e is 
minimized over the model parameter space. The results suggest that invariant minimization can 
correctly identify the model from the data. 



